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Abstract 

Rapidly growing numerical instabilities routinely occur in multidimensional 
particle-in-cell computer simulations of plasma-based particle accelerators, 
astrophysical phenomena, and relativistic charged particle beams. Reduc- 
ing instability growth to acceptable levels has necessitated higher resolution 
grids, high-order field solvers, current filtering, etc. except for certain ratios 
of the time step to the axial cell size, for which numerical growth rates and 
saturation levels are reduced substantially. This paper derives and solves the 
cold beam dispersion relation for numerical instabilities in multidimensional, 
relativistic, electromagnetic particle-in-cell programs employing either the 
standard or the Cole-Karkkainnen finite difference field solver on a staggered 
mesh and the common Esirkepov current-gathering algorithm. Good overall 
agreement is achieved with previously reported results of the WARP code. 
In particular, the existence of select time steps for which instabilities are 
minimized is explained. Additionally, an alternative field interpolation algo- 
rithm is proposed for which instabilities are almost completely eliminated for 
a particular time step in ultra-relativistic simulations. 

Keywords: Particle-in-cell, Esirkepov algorithm, Relativistic beam, 
Numerical stability. 
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1. Introduction 

In a laser plasma accelerator (LPA), a laser pulse is propagated through 
a plasma, creating a wake of very strong electric fields of alternating polarity 
PQ. An electron beam injected with the appropriate phase can thus be accel- 
erated to high energy in a distance much shorter than those for conventional 
acceleration techniques J2]. Simulation of a LPA stage from first principles 
using the Particle-In-Cell (PIC) technique in the laboratory frame is very de- 
manding computationally, as the evolution of micron-scale laser oscillations 
needs to be followed over millions of time steps as the laser pulse propagates 
through a meter-long plasma for a 10 GeV stage [3]. 

A method recently was demonstrated to speed up full PIC simulations 
of a certain class of relativistic interactions by performing the calculation 
in a Lorentz boosted frame [4], taking advantage of the properties of space- 
time contraction and dilation in special relativity to render space and time 
scales (which are separated by orders of magnitude in the laboratory frame) 
comparable in a Lorentz boosted frame, resulting in far fewer computer op- 
erations. In the laboratory frame the laser pulse is much shorter than the 
wake, whose wavelength is also much shorter than the acceleration distance 
(Ataser < X wa ke < ^ acceleration) ■ In a Lorentz boosted frame co-propagating 
with the laser at a speed near the speed of light, the laser is Lorentz ex- 
panded (by a factor (1 + Vf) 7/, where 7/ = (l — vj) , Vf is the velocity 
of the frame, normalized to the speed of light). The plasma (now moving 
opposite to the incoming laser at velocity —Vf) is Lorentz contracted by a 
factor 7/. In a boosted frame moving with the wake (i.e., 7/ ~ 'Jwake), 
the laser wavelength, the wake, and the acceleration length are comparable 
(\i ase r < X wake ~ X acce i eration ) , leading to far fewer time steps by a factor 
(1 + Vf) 7j, hence far fewer computer operations [U E]. 

A violent numerical instability, associated with the plasma back-streaming 
at relativistic velocity —Vf in the computational frame, limited early at- 
tempts to apply this method to speedups ranging between two and three 
orders of magnitude [31 El Ej- Control of the instability was obtained via 
the combination of: (i) the use of a tunable electromagnetic solver and an 
efficient wide-band digital filtering method [8], (ii) observation of the ben- 
efits of hyperbolic rotation of space-time on the laser spectrum in boosted 
frame simulations [9], and (iii) identification of a special time step at which 
the growth rate of the instability is greatly reduced [8]. The combination of 
these methods enabled the demonstration of speedups of over a million times 
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[9]. The instability is described in some detail in [S]. 

In this paper, the analysis first reported in [10J, which introduced the 
concept of numerical Cherenkov instabilities, is generalized and extended 
to two dimensions. (Extension to three dimensions follows readily from the 
analysis presented below, and the same conclusions apply) The new analysis 
recovers the salient features of the instability described in [5] , including the 
existence of the special time step. Growth rates are calculated for various 
cases of ultra-relativistic drifting plasmas and shown to match closely the 
growth rates obtained using the PIC code WARP [llj. Additionally, an 
alternative field interpolation algorithm is proposed for which instabilities 
are almost completely eliminated for a particular time step. A similar type 
of instability was reported in the calculation of astrophysical shocks [12], and 
the conclusions from this paper should apply readily. 

A general derivation of the numerical instability dispersion relation for 
multidimensional PIC codes employing the Esirkepov algorithm is outlined 
in Sec. 2. The dispersion relation is specialized in Sec. 3 to a cold, relativis- 
tic beam in two dimensions for comparison with WARP simulations. Sec. 
4 provides a simple yet reasonably accurate analytical expression for max- 
imum numerical instability growth rates and, thereby, identifies time steps 
for which growth rates are significantly reduced, or even eliminated under 
some conditions. Then, the dispersion relation is solved numerically for a 
range of parameters and compared with WARP results in Sec. 5. (Most 
of these analytical and numerical calculations were performed using Math- 
ematica [33]). Finally, Sec. 6 presents WARP simulations, demonstrating 
the near absence of numerical instabilities for an appropriately chosen field 
interpolation scheme and time step in two and three dimensions. 

2. Numerical instability dispersion relation 

The derivation here follows closely that of the general numerical insta- 
bility dispersion relation in [33], and only those steps that differ will be 
presented . To start, the present derivation is based on the electromagnetic 
fields themselves rather than on the potentials. 



dE 

~dt 



V x B - J, 



(1) 



dB 

~dt 



-V x E. 



(2) 
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Units are chosen such that, without loss of generality, the speed of light and 
other constants are unity. If the differential equations are replaced by cor- 
responding finite difference equations, and the difference equations Fourier- 
transformed in space and time, we obtain expressions of the form, 

HE = -[k]xB + iJ, (3) 

[u]B = [k] x E. (4) 

Brackets around quantities designate their finite difference representations. 

The Esirkepov algorithm determines not the current itself but its first 
derivative; see Eq. (19) of [15J. The Fourier transform of this equation can 
be written as, 

W x ] ( [k x ]J x ] 

W y \ = -iAt\ [k y ]J y \, (5) 
W z J ( [k z ]J z J 

with CT the current contribution of an individual particle, and At the simu- 
lation time step. W is further defined in terms of the current interpolation 
function S J by Eq. (23) of |15j . which when Fourier-transformed becomes, 



W x ) ( sin (k' x v x f) [cos (k' y v y f) cos (k' z v z f) - § sin (k' y v y f) sin (k' z v z f )] ] 

W y \ = -2iS J I sin (k' yVy f) [cos (k' z v z f) cos (k' x v x f) - | sin (k' z v z f) sin (k' x v x f)\ \ , 
W z J { sin (k' z v z f) [cos (k' x x y f) cos (k' y v y f ) - § sin (k' x v x f) sin (k' y v y f )] J 

(6) 

with v the particle velocity. Combining Eqs. ^ and ^ provides the desired 
expression for the particle current, 



J x ) ( sin (k' x v x f) [cos (k' y Vyf) cos {k' z v z f) - | sin (k' y v y f) sin (k' z v z f )] /[k x ] 

J y > = S J — { sin (k' y v y f) [cos ) cos (k' x v x f) - | sin (k' z v z f ) sin (^f )] /[^] 

j; J A * [ sin (fc> 2 f ) [cos (A;ix v f ) cos {k' y v y f ) - | sin (fc> x f ) sin (fe^f )] /[^] 

(7) 

This expression must, of course, be integrated over the linearized particle 
distribution function to obtain the total current. Note that Eq. ^ reduces 
to J7" = v in the limit of vanishing time step and cell size, as it should. 
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Modeling relativistic simulations requires replacing dv/dt by dp/dt, with 
p = 7v the relativistic momentum and 7 the relativistic energy, in Eq. (16) 
of [13] ■ This change flows through to Eqs. (17) and (18), in which 

— = - — - — p — (8) 
dp 7 <9v 7 3 <9v 

replaces d/dv. The force in Eqs. (17) - (19) of [TJ] applied to the particles 
is E + vxB, with the components of E and B multiplied by the Fourier 
transforms of their respective interpolations functions: 

S E *E x + v y S B *B z -v z S B *B y ) 
S E yE y + v z S B *B x - v x S B *B z \ . (9) 
S E 'E s +v a S B *B v -v v S B 'B x J 

In contrast to [H], the Fourier-transformed field interpolation functions are 
not assumed to be identical. 

Replacing appropriate parts of Eqs. (18) and (23) of [H] by the corre- 
sponding terms from Eqs. (jTj) - ([9]) yields 



\ F -\ 






hi 







£/ F 'l JCSC 



{u - k • v) — 



^/d 3 v (10) 



summed over spatial aliases m z and defined in [13]. The determinant 

of the 6x6 matrix comprised of Eqs. ([3]), Q, and (JTo]) is the desired dispersion 
relation. A striking difference between this and the general dispersion relation 
in [T3] is that the present dispersion relation contains trigonometric functions 
involving particle velocities. 



3. WARP 2-d dispersion relation 

For comparison with WARP two-dimensional, cold beam simulation re- 
sults |8j, we reduce Eqs. (3~j) and Q to a 3x3 system in {E z , E x , B y } and 
perform the integral in Eq( 10 ) for a cold beam propagating at velocity v in 



the ^-direction. The resulting matrix equation is 




£,z,x £>z,y ~h fox] 
ix,x + H £x,y ~ [K] 

-D* z [k z ] [u] 




(11) 
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D* and D* are introduced at this point to accommodate the Cole-Karkkainnen 
field solver, sometimes used in WARP; it is discussed near the end of this 
section. The quantities £ are employed purely for notational simplicity. 



E * esc 2 





At 




- k z v) — 


-sin 











K/%1 (12) 



^ = -n^2S J S 1 



CSC 





cot 






At / 




















2 ; 2 


T si H 


2 2 y 



k x /[k z 



(13) 



= fit; 5 J 5 B " esc 



cot 



At 

T 



sin ^/^J' 
(14) 



£x,x — n *y ] s j s 



CSC 







At / 






k'v) — 


— cos 


k'v — 




z j 2 


2 V 


z 2 J 



nv 



CSC 



k' z v] 



At 



Y cos \ k '* v 2 ) (U,) 



summed over spatial aliases, k' z = k z + 2ir/Az and fc^. = fc^ + m a . 2ir/Ax, 
with m 2 and m x integers. The resonances, u — k' z v, introduce an infinity of 
spurious beam modes with effective charge densities proportional to S J S Ez , 
etc. n is the beam charge density divided by 7, which can be normalized to 
unity. However, explicitly retaining it in the dispersion relation sometimes 
is informative. 

WARP employs the usual staggered spatial mesh and E-B leapfrog in 
time [16]. Hence, 

At\ f At 



\tu\ = sin \ us 



(17) 
(18) 



6 



sin ( k 



Ax 



(19) 



Also as usual, WARP employs splines for current and field interpolation. The 
Fourier transform of the current interpolation function is 



sin 



,Az 



I k 



, Az 



sm 



, Ax 



, Ax 



(20) 



and £ x are the orders of the current interpolation splines in the z- and 



rc-directions. So, for instance, an exponent of 2 in Eq. (20) corresponds to 



linear interpolation, and of 4 to cubic interpolation. Analogous definitions 
apply to the three field interpolation functions, but the spline orders need 
not be the same. WARP typically employs field interpolation splines like 
those of the currents but with the E z splines one order lower in z, the E x 
splines one order lower in x, and the B y splines one order lower in both. (This 
particular choice of spline orders is derivable by Galerkin's method |TT] and 
has superior energy conservation properties [HI [191 EQ]. It will be referred 
to subsequently as "Galerkin field interpolation".) 



S 1 



sin ( k 



, Az 



, Az 



(21) 



s 



•i) r 



(22) 



S 



COS u 



At 



sm 



k' 



Az 



/Ik 



Az 



sm 



Ax 



I k 



Ax 



-1) 



(23) 

The alias phase factors appearing at the ends of Eqs. (21 ) - (23) arise from 
the half-cell offsets from the current interpolation mesh of the corresponding 
fields. Averaging By in time before applying it to particles causes the factor 
cos (w^) in Eq.Q. 

Another credible choice of field interpolation functions is splines of the 
same order as those for the current interpolation function, in which case 
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Eqs. (21) 



(23) contain only powers of i x + 1 and £ z + 1. The powers of 
-1 are unchanged. This seemingly minor change has a significant impact on 
numerical stability for some choices of At. (It will be referred to subsequently 
as "uniform field interpolation".) 

The Cole-Karkkainnen field solver [2"T | |2"2" | |2"3"] . mentioned above, increases 
the Courant limit on the simulation time step and in some cases reduces 
numerical dispersion in the electromagnetic fields. It is discussed in some 
detail in Sec. 2.2 of [8J. For our purposes, 



D* 



A/3 X sin 2 



(24) 



D* = l- A/3 Z sin 2 k 



Az 



(25) 



For Ax = Az, the choice /3 X — /3 Z — 1/8 relaxes the Courant limit to At < 
Az, while minimizing numerical dispersion in the vacuum fields along major 
axes. 

Finally, we note that m x alias terms in the dispersion relation can be 
summed explicitly by means of Eqs. (1.421.3) and (1.422.3) of [21] or deriva- 
tives thereof, once choices have been made for the interpolation functions. 
For example, the fc^-dependent terms in £ zz sum to 



sin i K^pj i Oi^j 



2 cos ( k 



Ax 



/3 (26) 



for i x = 1. Note that the m x = term alone has the value (2/7r) 4 for k x 
near its maximum value, ir/ Ax. In contrast, the sum has the value 1 /z there. 
(Most of the difference is due to the m x = —1 alias, which is typical.) Since, 
as we shall see, peak growth rates typically scale as the cube root of such 
sums, the difference in predicted peak growth rates is of order 20%. 



4. Approximate peak growth rates 



£ ZiZ , defined in Eq. (12), scales as 7 (with n held constant) and can 
be ignored for highly relativistic calculations, on which this paper focuses. 
Likewise, 1 — v ~ 7~ 2 /2, and can be set to zero. Additionally, 



£,z,x£,x,y S2,2/S. 



(27) 
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is satisfied for individual modes and is satisfied approximately for cross- 
products between modes. With these assumptions the dispersion relation 
(the determinant of Eq. ([IT])) has the form 



Cq + n G\ esc 



, At 



£c 2 



n > ^ 2 esc 

in. 



with Co the vacuum dispersion function, 



C Q = [co] 2 -D* z [k x f-Dl [k z ]\ 



and 



Vx S J S E * cos I fc 



M 2 f [**] 



, At 



^ N 2 2^ k * b 



= 0. (28) 
(29) 

S E * S B y' 
(30) 



At 



c 2 = ^[fcJ T - £a£s 



5 i 



5 



cos 



At 

T 



sm 



At s 



Eq.(28) reduces, of course, to C 



n 



(31) 

in the limit of vanishing time step 
and cell size. All the beam modes in Eq.(28) are numerical artifacts, even 
the m z = mode. 

Coupling between these beam numerical modes and electromagnetic modes 
(the roots of Co = 0) gives rise to what has become known as the numerical 
Cherenkov instability [TU1 |2S], which can be quite virulent. Fig. [l] is a typi- 
cal normal mode diagram, showing the two electromagnetic modes and beam 
aliases m z = [-3, 3] for At = 0.7 Az, (3 X = (3 Z = 0, and k x = (Unless 
otherwise noted, other parameters for this and other figures are n = 1 and 
Ax = Az = 0.3868.) Fig. [2] depicts the locations in /c-space of normal mode 
intersections, such as those in Fig. [TJ as k x is varied. 

Comparing Fig. [2] with corresponding WARP results in Fig. [3] indicates 
that the strongest instabilities lie along the m z = -1 and resonance curves at 
larger k x . (The WARP simulations were performed on a 128 x 128 square grid 
with periodic boundary conditions and a uniformly distributed plasma mov- 
ing axially at an energy of 7 = 130, seeded with a small random transverse 
velocity. Plots similar to Fig. [3] appear in Also visible, although 

just barely, are much more slowly growing instabilities along the m z = +1 
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and m z = -2 resonance curves. We now proceed to estimate these instability 
growth rates. 

Resonance curves, such as those in Fig. [2j are given by Eq. (29) with to 
replaced by k' z , solved for k x as a function of k' z . (Recall that sm(k' z f) = 
sin 2 (^f ).) 



Ax 



arcsm 



A 



(f) 2 sin 2 (^f)-(|f)%in 2 (^) 



\ 



(32) 
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-4sin 2 (k' z f) (/3, + /3 2 (tf) 2 
To obtain an estimate of the numerical instability growth rate along a res- 



onance curve, we expand Co and the cosecants in Eq. (28) to first order in 



(ui — k' z ), set C\ = 0, and set u = k' z in The resulting cubic equation has 
one unstable root, 



Im (u) 



v5 
2 



\ 



2 u 



^S B v 



sin 



(Kf) 



sin 



(Kf) 



esc ( k 



evaluated at k T 



(33) 

Although it may appear that Eq. (33 ) becomes singular 
when k' z approaches zero, k r x approaches zero there also, as k' z 2 . Consequently, 
the growth rate vanishes in that limit. 

For completeness, we note that instability also occurs off-resonance when 
C0C2 > Ci /4, evaluated at u ~ k' z and arbitrary k x . The resulting growth 
rate is 



Im (lo) 



y/C?/4 - C C 2 

Co 



(34) 



Although off-resonance growth is weaker than on-resonance, it often occurs at 
smaller k z , where it may be more difficult to filter. (The residual instabilities 
after digital filtering discussed in the fourth paragraph of Sec. 5 are, for 
instance, of this sort.) 

Fig. [4] displays maximum instability growth rates for the Galerkin field 
interpolation algorithm as At/Az varies over its range of allowed values for 
0z = Ab (collectively, (3) = and Ys. (Intermediate values of (3 produce curves 
intermediate in shape.) The pronounced dip in both curves, at At/Az^i 0.66 
for (3 = and 0.69 for (3 = 1 /s, previously has been observed in simulations 
ED- It occurs because Im(uj) vanishes for some value of k 7 , which occurs 
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when 

At , { ,,Az\ ,,Az 



Az 



^^y) = fc^ sin (A£At) . (35) 



Eq. (35) has solutions for the m z = -1 and resonances only over a narrow 
range of time steps, a/2/2 > At/ Az > 0.65. Precisely where the minimum 
falls within this range depends on algorithmic details. 

Similarly, Fig. [5] displays maximum instability growth rates for the uni- 
form field interpolation algorithm as At/ Az varies over its range of allowed 
values for = and 1 /s. For all values of (3, the growth rate vanishes at 
At/ Az = 1 /2. Why this should be so is evident from 

^sin(^) =^sin(^At), (36) 

which differs from its Galerkin counterpart, Eq. (35 ), by a factor of sin (k' z ^f-) / (k' z ^f) 



Eq. (36) is satisfied for all k' z at At/ Az = 1/2, and for no values (apart from 



0) of k' z otherwise. 



Eq. (33) also provides a simple means for estimating the effect of current 
filtering on numerical Cherenkov instabilities, because the Fourier transform 
of the digital filtering function appears simply as a factor multiplying n. 
Given the substantial growth rates of this instability, filtering must reduce 
currents in regions of fc-space where the instability is strong by some three 
orders of magnitude. Of course, any physical phenomena occurring in those 
same regions also will be suppressed. Using higher order interpolation (i.e., 



larger E's in Eqs. (20) - (23)) also reduces numerical instability growth, 
especially for higher order aliases. However, for typical simulation parameters 
it reduces the m z = -1 and instability growth rates by comparable, modest 
factors. Employing cubic rather than linear splines, for instance, would be 
expected to reduce maximum growth rates by of order (2/n)^ 3 . On this 
basis current digital filtering usually is more cost effective than higher order 
interpolation for suppressing numerical Cherenkov instabilities. 



5. Numerical solutions 



Reliably finding the roots of Eq. (11) can be accomplished as follows. 
Given how strongly even linear interpolation suppresses all but the first few 
aliases, we safely can truncate the infinite series in m z to a range of, say, [-3, 
3]. (Indeed, the smaller range [-1, 0] works fairly well in most cases.) Then, if 
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the aliases are well separated in u — k space (as they are in, for instance, Fig. 
[T|, the growth rates for any particular alias can be evaluated with reasonable 
accuracy by expanding the dispersion relation as a fourth-order power series 
in (u — k' z v) for the k' z in question and calculating all roots with a polynomial 
root finder. On the other hand, if aliases are separated in frequency by only 
a few times the typical growth rates, the expansion converges slowly, and 
an iterative solution is required. The Mathematica p3] FindRoot routine 
was used for the results that follow in this section, with three real roots or 
one real root and one conjugate pair of roots found per alias. Evaluations 
were performed on a 65x65 array in fc-space, consistent with the 128x128 
spacial grid used in WARP for comparable simulations. Obtaining results 
for a typical set of parameters required about 15 minutes on a 2.8 GHz, 2 
processor desktop computer. 

Fig. [6] presents numerical growth rate predictions corresponding to the 
WARP results in Fig. [3] The m z = -1 alias dominates the growth spectrum 
with a maximum growth rate of 0.56 at short wavelengths in x. (The ap- 
proximate growth rate based on the analysis in the previous section is 0.48.) 
Also visible is the fast growing m z = alias. The much weaker m z = -2 and 
+ 1 aliases are evident at smaller k z . As noted in the previous section, the 
m z = -1, 0, and +1 aliases all can be seen in Fig. [3j although the last of these 
aliases is faint, consistent with its relatively slow growth. Fig. [7] depicts grow 
rates measured in this WARP simulation (actually the average of one hun- 
dred such simulations). The agreement between Figs. [6] and [7] is very good, 
especially when one considers the difficulty in measuring smaller growth rates 
in simulations, where nonlinear mode coupling and thermal noise can be sig- 
nificant. Thus, the method used to determine automatically the growth rates 
in WARP works well for the largest growth rates, which are of most interest 
in any particular simulation, but not so well for the smallest growth rates. 

Maximum numerical growth rates observed in WARP for the Galerkin 
and uniform current interpolation algorithms with /3=0 and 1 /s are com- 
pared with the predictions of linear theory in Figs. [8] and [9] Agreement 
between theory and simulation is very good. Qualitative agreement with the 
analytical estimates of the previous section is quite acceptable. The sudden 
rise of growth when At/Az nears unity for /3 = 1 /s comes from an instability 
of the field solver algorithm at the Nyquist limit and is mitigated by using 
one or more passes of bilinear filtering of the current density, as explained in 
Appendix A of [8] and shown below. 

Fig. [TU] illustrates the effects of digital filtering and of higher order inter- 
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polation, in this case ten passes of the bilinear filter (including two compen- 
sation steps) described in [8] , cubic or linear interpolation in z with Galerkin 
field gathering, and /3 = The digital filter has the effect of multiplying n 
in the dispersion relation by 

cos 16 ( K^] ( 5 - 4 cos 2 ( k K ^-X) cos 16 ( k T ^] ( 5 - 4 cos 2 ( t,.^] 



2 ) \ \ 2 ) ) \ 2 J V V 2 / 

(37) 

It effectively eliminates numerical instabilities for k z Az/ir > 0.2 or k x Ax/it > 
0.2. With linear interpolation the m z = alias dominates the numerical 
instabilty growth rate except in the vicinity of At/Az w 0.69, where the 
m z = +1 alias dominates. Growth rates are reduced by roughly a factor 
of four compared to those in Fig. |8j Cubic interpolation has negligible 
effect on the m z = alias but almost completely suppresses the m z = +1 
alias. The minimum growth rate, now at At/Az « 0.70, drops by a further 



factor of three. (Measuring the WARP instability growth rates for Fig. 10 
was particularly challenging due to competition between the weak numerical 
instabilities, and thermal and nonlinear effects.) 



As a further comparison between linear theory, Fig. 11, and WARP re 



suits, Fig. 12 (also averaged over 100 simulations), we present growth rates 
for Galerkin current interpolation with ^ = 0.69, and (3 = The dom- 
inant alias is m z = +1, occurring at the rather small axial wave numbers, 
1.5 < k z < 3.5, and at most k x values away from the fc z -axis. Modestly to 
the right is the m z = —2 alias, occurring at 3 < k z < 4 for large values of k x . 
Generally, we expect the m z = +1 and -2 aliases to have comparable growth 
rates, just as the m z = and -1 aliases typically do. Finally, the m z = —3 
alias is modestly above background on the far right. For all these modes, 
theory and simulation growth rates agree to within about 15%. However, a 
region of reduced growth rate in the band 5 < k z < 6 occurs only in Fig. 



12, although it can be produced in Fig. 11 by artificially removing the off- 
resonance m z = — 1 contribution. This minor discrepancy is apparent only 
for parameters very near those listed in this paragraph. 



6. Application to the modeling of laser plasma acceleration 

As a verification that the theory that has been developed in this pa- 
per applies to the modeling of LPAs, series of two and three dimensional 
simulations of a 100 MeV class LPA stage were performed, focusing on the 



13 



plasma wake formation, using the parameters given in table [TJ The velocity 
of the wake in the plasma corresponds to 7 ~ 13.2, and the simulations were 
performed in a boosted frame of 7/ = 13. 

Reference simulations were run in two and three dimensions for condi- 
tions where no instability developed, and the final total field energy Wf was 
recorded as a reference value in each case. Runs then were conducted for 
the Yee (/? = 0) and Cole-Karkkainnen (f3 = 1/8) solvers, with Galerkin and 
uniform field interpolations. The final energy Wf was recorded and divided 
by the reference energy W/o. The ratio Wf/W/a is plotted versus time step 
from two dimensional simulations in Fig. [13] and from three dimensional 



simulations in Fig. 14, using linear current deposition and no smoothing of 
current and fields. Following the theoretical predictions, for the Galerkin 
interpolation scheme the instability is minimal around At/Az « 0.65 when 
(3 = and around At/Az w 0.69 when /3 = 1/8, while for the uniform 
interpolation scheme the instability is minimal around At/Az « 0.5. The 
ratio Wf/Wfo also is plotted versus time step from two dimensional simula- 



tions in Fig. [15] and from three dimensional simulations in Fig. 16, using, 
as is common practice in the modeling of laser plasma stages, cubic current 
deposition and 1 pass of bilinear smoothing plus compensation of current 
and fields gathered onto macroparticles. The beneficial impact on stability 
of smoothing and high order deposition is evident from the relatively wide 
band of stability that is available around At/Az ~ 0.5 with uniform gather, 
and the narrower band of stability that is available around At/Az m 0.7 
with Galerkin gather. This verifies that the theoretical results apply to real 
case simulations in two and three dimensions. 



7. Conclusion 

The numerical stability properties of multidimensional PIC codes employ- 
ing the Esirkepov current algorithm have been derived. Just as in PIC codes 
employing earlier current algorithms, here also fast-growing numerical insta- 
bilities are predicted for relativistic beam simulations. These instabilities 
can, of course, be reduced significantly by short wavelength digital filter- 
ing. However, time steps have been identified at which instability growth 
is reduced even without filtering. Particularly noteworthy is uniform field 
interpolation with At/Az = 1 /2 and any value of j3, for which simulations are 
numerically stable in the large 7 limit. These results have been confirmed 
with the WARP simulation code. 
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Table 1: List of parameters for a LPA stage simulation at 100 MeV 



plasma density on axis 


n e 


10 19 cm 3 


plasma longitudinal profile 




flat 


plasma length 


L p 


1.5 mm 


plasma entrance ramp profile 




half sine 


plasma entrance ramp length 




20 fim 


laser profile 




do exp (— r 2 /2a 2 ) sin (tiz/3L) 


normalized vector potential 


a 


1 


laser wavelength 


A 


0.8 /im 


laser spot size (RMS) 


a 


8.91 /im 


laser length (HWHM) 


L 


3.36 /im 


normalized laser spot size 


kp<j 


5.3 


normalized laser length 


kpL 


2 


cell size in x 


Ax 


A/32 


cell size in y (3D only) 


Ay 


A/32 


cell size in z 


Az 


A/32 


# of plasma particles/cell 




1 macro-e~+l macro-p + 



Additionally, WARP LPA simulations performed using uniform field in- 
terpolation with At/ Az = l /2 have demonstrated the practical value of this 
choice of parameters in two and three dimensions. The uniform field inter- 
polation offers much reduced growth rates, enabling faster simulations with 
fewer grid cells, lower order interpolation, and reduced digital filtering. In 
three dimensions, it enables existing PIC codes that incorporate the Yee 
solver, but not the CK solver, to benefit from the reduced growth rates at 
the special time steps over a wider range of cell aspect ratios (for cubic cells 
for example, the special time step is accessible only to the CK solver for 
Galerkin gather, while it is accessible to both the Yee and the CK solvers 
for uniform gather). The results that were obtained here also should apply 
readily to more efficient modeling of astrophysical shocks that use the same 
algorithms. 

Finally, the salutary effect of trigonometric functions involving particle 
velocities in the dispersion relation of the Esirkepov algorithm suggest that 
further improvements in PIC code stability can be achieved by developing 
field interpolation algorithms that introduce similar trigonometric functions, 
perhaps along the lines of Sec. 4 in [T4"] . 
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Figure 1: Normal mode diagram for^| = 0.7, /3 = 0, and k x — 1 /2 showing numerically 
distorted electromagnetic modes and spurious beam modes, m z = [—3, 3]. Numerical 
Cherenkov instabilities occur near mode intersections. 
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Figure 2: Locations in fc-space of resonances between electromagnetic modes and beam 
modes, m z — [—3, 3] for ^rf = 0.7 and /3 = 0. Intersecting resonance curves occur at 
different frequencies and, therefore, do not interact. 
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kz 

Figure 3: Fourier-transformed E z (log scale) at t = 16 from a WARP simulation with 
Galerkin field interpolation, ^ = 0.7, and (3 = 0. 
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Galerkin Field Interpolation 
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Figure 8: Maximum numerical instability growth rates observed in WARP and calculated 
from the numerical dispersion relation for Galerkin field interpolation with f3 = 0, 
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Uniform Field Interpolation 




At/Az 



Figure 9: Maximum numerical instability growth rates observed in WARP and calculated 
from the numerical dispersion relation for uniform field interpolation with f3 = 0, 
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Galerkin p=Ya Interpolation Filtered 
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Figure 10: Maximum numerical instability growth rates observed in WARP and calculated 
from the numerical dispersion relation for digital filtering as described in Sec. 5, overall 
linear or cubic interpolation in z, and Galerkin field interpolation with (3 = 1 /r. 
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Figure 13: Field energy relative to stable reference level vs At/Az from two dimensional 
WARP LPA simulations at 7 = 13, using Galerkin and uniform field interpolation with 
(3 — 0, 1 /s 1 no filtering, and linear interpolation. 
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Figure 14: Field energy relative to stable reference level vs At/Az from three dimensional 
WARP LPA simulations at 7 = 13, using Galerkin and uniform field interpolation with 
P = 0, Vs. 
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Figure 15: Field energy relative to stable reference level vs At/ Az from three dimensional 
WARP LPA simulations at 7 = 13, using Galerkin and uniform field interpolation with 
P = 0, Vs. 
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Figure 16: Field energy relative to stable reference level vs At/Az from three dimensional 
WARP LPA simulations at 7 = 13, using Galerkin and uniform field interpolation with 
= 0, Vs. 
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